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Abstract 

The pseudo-polar Fourier transform is a specialized non-equally spaced Fourier transform, 
which evaluates the Fourier transform on a near-polar grid known as the pseudo-polar grid. 

The advantage of the pseudo-polar grid over other non-uniform sampling geometries is that 
the transformation, which samples the Fourier transform on the pseudo-polar grid, can be 
inverted using a fast and stable algorithm. For other sampling geometries, even if the non- 
equally spaced Fourier transform can be inverted, the only known algorithms are iterative. 

The convergence speed of these algorithms as well as their accuracy are difficult to control, 
as they depend both on the sampling geometry as well as on the unknown reconstructed 
object. In this paper, a direct inversion algorithm for the three-dimensional pseudo-polar 
Fourier transform is presented. The algorithm is based only on one-dimensional resampling 
operations, and is shown to be signihcantly faster than existing iterative inversion algorithms. 

Keywords: 3D pseudo-polar Fourier transform, Radon transform, Unequally spaced FFT, Polar 

Fourier transform, Toeplitz matrices 

AMS subject classification: 65T50, 44A12, 92C55 

1 Introduction 

The fast Fourier transform (FFT) is one of the most commonly used algorithms, and has far 
reaching implications in science and technology [7j. While the FFT efficiently computes the discrete 
Fourier transform (DFT) on a Cartesian grid, in many applications it is required to compute the 
Fourier transform on non-Cartesian grids. Computing the FFT on a non-Cartesian grid can be 
implemented using one of the available non-equally spaced FFTs (NUFFTs) [HI HDl HSl HH l36[ 120]. 
Of particular importance in applications are the polar grids in two and three dimensions, which 
emerge naturally in a wide range of applications, from computerized tomography to nano-materials 
to image processing [aHiEiiiiiiiiiEsiianiEHiETiiiiEn]. 
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Although the Fourier transform can be evaluated on a polar grid by using one of the available 
NUFFT algorithms, there are two factors that may limit their applicability to large problems. 
First, although all NUFFT algorithms have the same asymptotic complexity as the classical FFT, 
namely, 0{n\ogn), where n is the number of input/output samples, their runtime complexity 
involves rather larger constants. This is due to local interpolations used by all these methods 
to resample from the Cartesian grid to the polar grid. Therefore, Althongh these algorithms are 
asymptotically as efficient as the FFT, they are prohibitively slow for large input sizes. Second, 
as many real-life problems can be formulated as the recovery of image samples from frequency 
samples (e.g., medical imagery reconstrnction algorithms), they actnally reqnire compnting the 
inverse Fonrier transform. However, unlike the FFT, whose inverse is eqnal to its adjoint and thus 
can be easily inverted, the NUFFT, except for very special cases, does not have this property, and 
therefore, more elaborated inversion schemes are reqnired. In particular, although the transform 
that evalnates the Fourier transform on the polar grid is formally invertible, the condition nnmber 
of this transformation exclndes snch inversion in practice. 

Although the polar sampling geometry of Fourier space is natural for various applications snch 
as compnterized tomography [9] or electron microscopy [MIE], h is sometimes possible to overcome 
the aforementioned difficulties by replacing the polar freqnency sampling grid with the psendo- 
polar freqnency sampling grid. The psendo-polar grid [2] (to be described below) is nearly polar, 
but unlike the polar grid, it is not equally spaced in the angular direction, but uses equally spaced 
slopes. We refer to the transformation which samples the Fourier transform on the psendo-polar 
grid, as the psendo-polar Fonrier transform (PPFT). Both the PPFT and its inverse admit efficient 
and stable algorithms. Iterative algorithms nsed to invert the psendo-polar Fourier transforms can 
be fonnd in [21EU [121 [37] . Since the psendo-polar Fonrier transform is ill-conditioned, these inver¬ 
sion algorithms reqnire preconditioning. Several preconditioning approaches have been suggested, 
such as using Voronoi weights [13] or a preconditioner derived from the Jacobian [3^. In partic- 
nlar, it is shown in [39] that nnder the appropriate preconditioning, the condition nnmber of the 
pseudo-polar Fourier transform is small. Therefore, it can be efficiently inverted using iterative 
inversion schemes. Nevertheless, snch an iterative inversion reqnires to compute the pseudo-polar 
Fourier transform and its adjoint in each iteration, which as mentioned above, is too expensive for 
large problems. Even if the psendo-polar Fourier transform and its adjoint are computed using 
specialized algorithms [3], the number of iterations reqnired for the inversion may be of the order 
of a few dozens, which may be too slow for large scale problems. Moreover, the accuracy of the 
reconstructed object by such iterative methods typically cannot be estimated from the convergence 
criterion and more calculations are needed to verify that the reconstructed object has the required 
accuracy. 

The PPFT can also be expressed as a multilevel Toeplitz operator EH. Thus, the inverse PPFT 
is eqnivalent to the inversion of this operator. Direct and iterative algorithms have been devised 
for inverting Toeplitz-related operators. Direct inversion algorithms, which are stable and fast. 
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are available for Toeplitz matrices inversion mn- On the other hand, these algorithms do not 
extend to multilevel-Toeplitz operators. The only available exact algorithm for inverting two-level 
Toeplitz operators |13] is not applicable to our setting since it assumes that the two-level Toeplitz 
matrices are triangular in one or both their levels. 

In this paper, we present a direct algorithm for inverting the three-dimensional pseudo-polar 
Fourier transform, which generalizes the two-dimensional direct inversion algorithm in [2]. Unlike 
the iterative approaches, the proposed algorithm terminates within a fixed number of operations, 
which is determined only by the size of the input. The inversion algorithm is shown numerically to 
have high accuracy, and since most of its steps are implemented in-place, its memory requirements 
are sufficiently low to be applicable to very large inputs. Finally, the algorithm is numerically 
stable since it consists of a series of well-conditioned steps. In particular, no preconditioner is 
needed. 

This paper is organized as follows. In Section we provide the required background on 
the pseudo-polar grid and the pseudo-polar Fourier transform. In Section we give a detailed 
description of the direct inversion algorithm for the three-dimensional pseudo-polar transform. 
Numerical results and comparisons to other algorithms are given in Section Finally, some 
concluding remarks are given in Section 


2 Mathematical preliminaries 


In this section, we provide the mathematical background required to derive our inversion algorithm. 
In Section [2.1| we revisit the pseudo-polar Fourier transform. In Section 2^, we describe a fast 


algorithm for solving Toepliz systems of equations. This algorithm is used in Section 2^ for fast 
resampling of univariate trigonometric polynomials. 


2.1 Pseudo-polar Fourier transform 

The 3D pseudo-polar grid, denoted flpp, is defined by 

Dpp ^ppx D ^ppy U ^ppz , 


( 1 ) 


where 


'ppx 


•ppy 


•ppz 






( 2 ) 


with m = qn + lioi an even n and a positive integer q. We denote a specific point in ^ppx, ^ppy and 
^ppz by flppx{k, l,j), flppy{k, l,j) and Qppz{k, l,j), respectively. The psuedo-polar grid is illustrated 
in Fig. for n = 4 and q = 1 . As can be seen from Fig. the pseudo-polar grid consists of 
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Figure 1: The three sectors Vtpp^^ ^ppx and Q.ppz that form the 3D pseudo-polar grid. 


equally spaced samples along rays, where different rays have equally spaced slopes but the angles 
between adjacent rays are not equal. This is the key difference between the pseudo-polar grid and 
the polar grid [1] . Thus, we can refer to fc as a “pseudo-radius” and to I and j as “pseudo-angles”. 

Next, we dehne the discrete time Fourier transform of an n x n x n volume I by 

nl2-l 

G [-m/2,m/2], (3) 

u,v,w=—nl2 

where as before m = qn + 1 for an even n and a positive integer q. The three-dimensional 
pseudo-polar Fourier transform (PPFT) is defined as the samples of the discrete time Fourier 
transform of Eq. on the pseudo-polar grid Qpp. Specifically, if we denote the concatenation of 
three arrays Ai, A 2 , and A^ by A = [Ai, A 2 , A^], then the pseudo-polar Fourier transform of a 
volume I G is an array Ifipp G C^xD+i)xD+i)x(3«+i) given by 

In 

i Lpp 


1 Iq 1 Id 

z ; i C‘ppy 7 i ' 


(4) 


where 


= n k,--k,- — k] , 

^‘‘vvx I 5 Tl ^ Tl ^ 








and Dppa,, and flpp^ are dehned in Eq. The parameter q in Eqs. and determines the 
frequency resolution of the transform, as well as its geometric properties. For example, as shown 
in [3], in order to derive a three-dimensional discrete Radon transform based on the pseudo-polar 
Fourier transform, q must satishes q > 3. 

The pseudo-polar grid appeared in the literature several times under different names. It was 
originally introduced by |35] under the name “Concentric Squared Grid” in the context of com¬ 
puterized tomography. More recent works in the context of computerized tomography, which take 
advantage of the favorable numerical and computational properties of the grid include [33l IM] , 
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where equally sloped tomography is used for radiation dose reduction. Other image processing 
applications that use the pseudo-polar grid include synthetic aperture radar (SAR) imaging [2H], 
Shearlets [271 EH], registration [32|, and denoising HU, to name a few. 

Recently, [39] proposed fast iterative inversion algorithms for the 2D and 3D pseudo-polar 
Fourier transforms, which are based on the well-known convolution structure of their Gram opera¬ 
tors, combined with a preconditioned conjugate gradients [T8| solver. The convolution structure of 
the transform allows to invert it using the highly optimized FFT algorithm instead of the forward 
and adjoint transforms derived in [2]. As the conjugate gradients iterations require precondition¬ 
ing, [39] proposes a preconditioner that leads to a very small condition number, thus making the 
inversion process fast and accurate. However, this approach has drawbacks such as high memory 
requirements and dependency of the number of iterations and on the size of the data. The algo¬ 
rithm presented in this paper, which has a low memory requirements, achieves high accuracy and 
is faster than the iterative algorithm [39] . 


2.2 Solving Toeplitz systems 

Let A„ be an n X n Toeplitz matrix and let y be an arbitrary vector of length n. We describe a fast 
algorithm for computing A~^y. This algorithm is well-known and appears, for example, in [2], but 
we repeat it here for completeness of the description. The algorithm consists of a fast factorization 
of the inverse Toeplitz matrix followed by a fast algorithm that applies the inverse matrix to a 
vector mi El]. We denote by T„(c, r) an n x n Toeplitz matrix whose first column and row are c 
and r, respectively. For symmetric matrices, c = r. 

Circulant matrices are diagonalized by the Fourier matrix. Hence, a circulant matrix Cn can be 
written as Cn = W*DnWn, where Dn is a diagonal matrix containing the eigenvalues Ai,..., A„ of 
Cn, and Wn is the Fourier matrix given by Wn{j, k) = Wg27njA:/n_ Moreover, if c = [cq, ci,..., Cn-i]'^ 
is the first column of Cn, then WnC = [Ai,..., A„]^. Obviously, the matrices Wn and W* can be 
applied in 0{nlogn) operations using the FFT. Thus, the multiplication of Cn with an arbitrary 
vector X of length n can be implemented in 0{n log n) operations by applying FFT to x, multiplying 
the result by Dn, and then applying the inverse FFT. 

To compute AnX for an arbitrary Toeplitz matrix An = Tn{c,r) and an arbitrary vector x, we 
first embed An in a 2n x 2n circulant matrix C 2 n 


Con = 


An Bn \ 
Bn An / 


where is an n x n Toeplitz matrix given by 


Bn = TndO, Tn-l, • • • , Ta, Ti], [O, C„_i, . . . , Ca, Ci]). 
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Then, AnX is computed in 0 (nlog? 7 ,) operations by zero padding x to length 2n, applying C 2 n to 
the padded vector, and discarding the last n elements of the resulting vector. 

Next, assume that An is invertible. The Gohberg-Semencul formula [El US] provides an explicit 
representation of A~^ as 

A-^ = - (M1M2 - M3M4), (6) 

Xo 

where 


Ml = T„([a;o,a;i,.. .,Xn-i], [a:o, 0 ,... , 0 ]), ( 7 ) 

M 2 = T„([?/„_i,0, ...,0],[|/„_i,|/„_2,...,|/o]), (8) 

M 3 = T„([ 0 ,?/o,---,l/n- 2 ],[ 0 , ..., 0 ]), ( 9 ) 

Mi = Tn{[0, . . . , 0 ], [ 0 , Xn-l, ..., Xi]), . ( 10 ) 


The vectors x = [xq, ..., Xn_i] and y = [|/o,..., l/n-i] are given as the solutions to 

AnX = eo, Co = [1,0,... ,0]'^, 

AnV 671—1) Sn—1 [0, • • • ) 0, 1] . 


The matrices Mi, M 2 , M 3 , and M 4 have Toeplitz structure and are represented implicitly using 
the vectors x and y. Hence, the total storage required to store Mi, M 2 , M 3 , and M 4 is that 
of 2n numbers. If the matrix is fixed, then the vectors x and y can be precomputed. Once 

the triangular Toeplitz matrices Mi, M 2 , M 3 and M 4 have been computed, the application of 
A~^ is reduced to the application of four Toeplitz matrices. Thus, the application of A~^ to a 
vector requires O(nlogn) operations. The pseudo-code of applying A~^ to a vector is described 
in Algorithms]^ and in the appendix. Algorithm lists the function ToeplitzDiag, which 
computes the diagonal form of the circulant embedding of a Toeplitz matrix. Algorithm lists 
the function ToeplitzInvMul, which efficiently multiplies an inverse Toeplitz matrix, given in its 
diagonal form, by a vector. The latter function uses ToeplitzMul, listed in Algorithm which 
efficiently multiplies a general Toeplitz matrix by a vector. 


2.3 Resampling trigonometric polynomials 

The main (and in fact the only) tool behind our algorithm is resampling of univariate trigonometric 
polynomials. Assume we are given a set of points yj G [—vr, tt], and their values {/(|/j)}^ 4 , 

where / is some unknown univariate trigonometric polynomial of degree n < N. We want to 
estimate the values of / at a new set of points {xj}^i, Xj G [—7r,7r]. This can be formulated by 
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first solving 


mm 

OL 


N 

i=i 


n/2-l 

/(»)- E 

k=—n/2 


a — (a_„/ 2 ,..., Q;n/ 2 -i) £ C"", 


( 12 ) 


for the coefficients vector a, followed by evalnating 


n/2-l 

/(i.)= E 

k=—nl2 


(13) 


In matrix notation, Eq. 12 is written 


min 11/- Aalla, (14) 

a 

where the entries of the matrix A are given by ajk = and the coordinates of the vector / are 
fj = fiVj) (we denote by / both the vector and the function, as the appropriate meaning is clear 
from the context). Direct solution for the coefficients vector a requires 0{n^) operations (assuming 
N and M are of size 0{n)). Obviously, solving for a also depends on the condition number of A, 
which affects the accuracy of f{xi). 

We next present a fast algorithm for computing f{xi) (faster than directly solving hrst for a), 
which exploits the Toeplitz structure of A*A and uses the non-equally spaced EFT (NUFFT) [HI 
oroi HSl [la [Ml EH]. The resulting algorithm has complexity of 0{n\ogn + nlog 1/e) operations, 
where e is the accuracy of the computations, in addition to a preprocessing step that takes 0{n‘^) 
operations. 

Following [in], we dehne the NUFFT type-I by 

rij 

/fc = — k = -n/2,...,n/2-l. (15) 

rii 

^ j=i 

and the NUFFT type-II by 




njl-l 

E 

k=—nl2 


zkikxj 


j = l,...,n,-. 


(16) 


Both types can be approximated to a relative accuracy e in 0{nlogn + n log 1/e) operations by 
any of the aforementioned NUFFT algorithms. From the dehnitions of NUFFT of types I and II, 


we see that A* f, where A is the matrix from Eq. 14 and / is an arbitrary vector, is equal to the 
application of NUFFT type-I to /. Similarly, the application of A to an arbitrary vector c is equal 
to the application of NUFFT type-II to c. Hence, the application of A or A* to a vector can be 
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implemented in 0{n\ogn + nlog 1/e) operations. 


To solve the least-squares problem of Eq. 14 we form the normal equations 


A*Aa = A*f. 


( 17 ) 


The right-hand side in Eq. 17 can be computed efficiently using NUFFT type-I. The matrix ^4*^4 
is a symmetric Toeplitz matrix of size nxn. The hrst column of A* A, which due to the symmetric 
Toeplitz structure encodes all its entries, is computed efficiently by applying NUFFT type-I to 
the vector w whose entries are Wj = _ Computing (^4*^4)“^ takes 0{n‘^) operations using 

the Durbin-Levinson algorithm |3T] and applying it to a vector takes 0{n\ogn) operations using 
the Gohberg-Semencul formula [23], as was described in Section 2.2 This procedure is described 
in details in [2]. Since {A*A)~^ depends only on n, it can be precomputed. Therefore, solving 


for a in Eq. 17 takes (!2(?7,logn-|-nlog 1/e) operations and computing fix^) in Eq. 13 takes addi¬ 


tional 0{n\ogn + nlog 1/e) operations using NUFFT type-II. The entire resampling algorithm is 
described in Algorithm 


3 Direct inversion of the 3D PPFT 

Given the PPFT dehned in Eq. of an unknown n x n x n volume /, the proposed direct 
inversion algorithm recovers I in two steps. The first step resamples the PPFT into an intermediate 
Cartesian grid. Specifically, we resample the trigonometric polynomial / of Eq. from the grid 
VLpp in Eq. to the frequency grid 

hlc = Qv, qw) : u,v,w = —n/2,..., n/2}. (18) 


The second step of our algorithm recovers I from the samples of I on flc. Note that this second 
step cannot be directly implemented by inverse FFT. The two steps of the algorithm are described 
in Sections]^ and|^ respectively, its pseudo-code is given in Algorithm]^ and its complexity is 


analyzed in Section 3.3 


3.1 Resampling the pseudo-polar grid to the grid Qc 

We start by describing the procedure for resampling / from flpp to tic- If is based on an “onion 
peeling” approach, which resamples at iteration k, k = n/2,... ,0, points in Qpp with “pseudo¬ 
radius” k to points in tic that lie on a plane. When the iteration that corresponds to /c = 0 is 
completed, the values of / have been computed on all points of the grid lie- We define the sets 
C fig, k = n/2 ,..., 0, consisting of the frequencies on which I has been resampled until (and 
including) iteration k. Note that for convenience, we index the iterations from k = n/2 to fc = 0. 
Thus, we have that C and = fig. Formally, we define 0*-^^ to be the frequencies on 














which we resample I at iteration k so that 


^(fe) _ (](fc+i) u0{fc)^ 


where 0^^^ is give by 

0(U = 0U) u 0(_U y 0 

U U 0f] u 0L"], 

(19) 

and for j, 1 = —k,... ,k 

0 +i ={{(lk,qj,ql)}, 

Q-l={{-qk,qj,ql)}, 



©S ={id,qk,qi)}, 

1 

© 

(20) 



(1? 

1 

© 



For k = n/2 we set Moreover, the frequencies of are the frequencies 

Qpp^{±n/2,j,l), Qppy{±n/2,j,l), and Qpp^{±n/2, j,l), l,j = -n/2,.. .,n/2, defined in Eq.|^ Thus, 
the values of J on are given as a subset of the values of J on flpp. Therefore, no resampling 

is needed at this iteration. 

For subsequent iterations, we use an example to accompany and clarify the formal description. 
We use as an example a grid corresponding to n = 8 that demonstrates how the values of I on 
0(n/2-i) jg 0(3) particular example) are recovered. Since the same procedure recovers 

the values of / on &±l, 0±y, and Q±l, we only explain the process for the grid 0+]. The pseudo¬ 
code for recovering the values of J on 0^^^ is given in Algorithm]^ whose details are also explained 
below. 

green solid squares correspond to points from on which J has already been 


In Fig. 


2 a 


evaluated (in previous iterations), red circles correspond to points from Qppx with a fixed k = 
n/2 —1, and blue dots correspond to points of 0+,^ on which we want to evaluate /. The frequencies 
0^,^ are points in however, as can be seen from Eq. 20, they all lie on the same plane, and 


m 


therefore, we depict them as a two-dimensional image whose axes are i and j from Eq. 

consists of resampling the points of (solid green 


The first step, depicted in Fig. 


2 b 


squares) to the same spacing as the pseudo-polar points. In the case of /c = n/2 — 1, this means 
that the first and the last rows are resampled, and the result of this resampling is denoted by 
patterned yellow squares. This step is implemented by lines 11 — 14 in Algorithm Next, as 
depicted in Fig. 3a, for all columns, we use the latter resampled points (patterned yellow squares) 
together with the points of the pseudo-polar grid at the same column to resample / to intermediate 


sampling points. Figures 3a and pm depict the resampling of one such column. Note that the 


resampled values of I after this step, which are depicted as hlled teal circles in Fig. are not yet 
the values of I on the points of Q^x\ Up to here, we only applied resampling to the columns (but 
not yet to rows). We repeat this for all the columns (see Fig. 4a for resampling of another column), 
which results in the grid of Fig. Tb This is implemented by lines 18 — 20 of Algorithm which 
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Figure 2: The grid before and after the applications of lines 11 — 14 in Algorithm]^ Green solid 
squares correspond to points from red circles correspond to points from ^ppx with a hxed 

(k) 

k = n/2 — 1, blue dots correspond to points of ©Y*, and patterned yellow points correspond to 
point of resampled to their intermediate position. 


using the conventions of Fig. resamples the red circles along the columns to the filled teal circle, 
which are now on the same rows as the points indicated by the blue dots (our target sampling 
points). Next, we apply one-dimensional resampling to each row, by using the resampled points 
from the previous steps together with the points of (Fig- 5a) to recover the samples of ©i^^ 

(Fig. 5b). This is implemented by lines 22 — 24 of Algorithm]^ All one-dimensional resampling 
operations in Algorithm]^ are implemented using the function ToeplitzResample in Algorithm]^ 


3.2 Recovering / from I on 

The second step in Algorithm (line 11) recovers the volume / from the samples in J on Gc 


(see Eq. 18). Define the operator Fd : C” —)■ C"'+^, by 


n/2-l 

{FDu)k = ^ k = -n/2,, n/2, m = qn + 1. (21) 

j=-n/2 

For a volume / of size n x n x n, we denote by Fj^\ F^^^ and F^'^ the application of Fn to the x, 
y and ^ directions of /, respectively. Furthermore, we dehne 


J — piv) pG) T 

1 — r jy r r P) i . 


( 22 ) 
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(a) Before resampling 


(b) After resampling 


Figure 3: The grid before and after the applications one iteration of lines 18 — 20 in Algorithm 



Figure 4: The grid before and after further iterations of lines 18 — 20 in Algorithm 
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Figure 5: The grid before and after the application of lines 22 — 24 in Algorithmj^ (row resampling). 


Algorithm 1 Direct inversion of the 3D pseudo-polar Fourier transform. 


1 : Input: Samples of the 3D pseudo-polar Fourier transform In^py, given by Eq. 

each of size {qn -|- 1) x (n -f 1) x (n -T 1). 

2: Output: Volume I of size n x n x n. 

Id ^ zeros(n -|- 1, u -|- 1, u -f-1) 

for j = 1 ,..., n/2 do 

{dU,-,-) <- recover2d(Jopp^(g(j -1) + 1,^(j,j - 1) > Algorithm 

6 : Inin + 2 - j,^ recover 2 d(Jnpp^(gn - q{j - 1) + l,n + 1 : -1 : l,n + 1 : -1 : 

l)Jp{n + 2--1) ^ 

7: ^ recover2d(Jopp^(g(j - 1) + 1,^(i, j,:), j - 1) 

8: /£,(:, u -1- 2 - j,:) ^ recover2d(4pp^(gn - g(j - 1) -f- 1, n + 1 : -1 : 1, n-f 1 : -1 : 

9: ^ recover2d(Jnpp^(g(j - 1) + 1,^(i,j), j - 1) 

10: Id{-, :,n + 2 - j) ^ recover2d(Jnpp,(gn - g(j - 1) -M, n-M : -1 : 1, n-t-1 : -1 : ■ 

,n + 2- j),j -1) 

11 : end for 

12: /£»(n/2 -f l,n/2 -f l,n/2 -f 1) ^ iQ^^^{3n/2 + l,n/2 + l,n/2 -f 1) (center point) 

13: / InvDecimatedFreq(//)) > (recover / from Id by Algorithm]^ 
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Algorithm 2 recover 2 d. This function recovers a 2D slice of the 3D pseudo-polar grid and 
recovers it to Cartesian grid (resample / from the set to the set 


1 : Input: - Slice of the 3D pseudo-polar Fourier transform of size (n -|- 1) x (u -|- 1). 

Ud - Already resampled slice oi Id from Algorithm of size (n -|- 1 ) x (n -|- 1 ). 
j - An integer indicating the stage within the peeling procedure. 

2 : Output: result of size (n -f-1) x (n - 1 - 1) (resampled part, Id from Algorithm]^ 

3: a ^ {nl2- j)l{nl2) 

4: m ^ qn + 1 
5: if j == 0 then 
6 : result <(— Un 

i Lpp 

return 
7: end if 

8 : yi -(r- [—n/2 : n/ 2 ] x (— 2 ) x q x Tr/m 
9 : Xi -ir- [—n/2 : n/2] x (—2) x q x a x tt/ m 
10 : C ^ zeros(n -|- 1, n -|- 1) 

11: for fc = 1,..., J do 

12: C(fc,:) ^ ToeplitzResample (17npp(fc, 1 : n-|-1), 1 / 1 , Xi j > Algorithm 5 


13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 


C(n -|- 2 — A:,:) ^ ToeplitzResample [U^ppin -|- 2 — fc, 1 : n -f l)?/i, Xi j 

end for 

X ■(— [—(n/2 — j) : {n/2 — j)] x (—2) x q x Ti/m 

y ^ [[~'^/2 : —n/2 + j — 1], [—n /2 : n/ 2 ] x a, [n /2 — j + 1 : n/ 2 ]] x (— 2 ) x g x vr/m 
Ri ■(— zeros(n -|- 1 — 2 j, n -|- 1 ) 
for A: = 1,..., n -f 1 do 

Ri ^ ToeplitzResample 
end for 

i ?2 ^ zeros(n -|- 1 — 2 j, n -R 1 — 2j) 
for A: = 1,..., n -|- 1 — 2j do 


C{1 : j, k), tJnppi-, k), C{n - j + 2 :n + l,k) 




i?2 t— ToeplitzResample 


UD{k + j,l : j),Ri{k,:),UD{k + j,n-j + 2 : n-M) 




end for 

result <(— JJd 

result{l + i : n - 1 - 1 — j, 1 -f j : n - 1 - 1 — j) ■(— i ?2 
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From the definitions of Fd, I and f2c in Eqs. and f^, respectively, we conclnde that I is 

eqnal to the samples of I on the grid ffc- Thus, if we apply the inverse of Fd to the x, y, and 2 ; 
dimensions of I (in a separable way), we recover I from the samples of I on f2c. The inversion of 
the operator Fd is described in [ 2 ]. Specihcally, in our case, for m = qn + 1, the adjoint of Fd is 
given by 

nl2 

(F», = j = -n/2,...,n/2-l. 

k=—nl2 

The matrix F^Fjj is a Toeplitz matrix, whose entries are given by 


n/2 

{F*r,FD)k,i= k,l 

j=-nl2 


—n/ 2 ,..., n /2 — 1 , 


and its inverse is applied as described in Eq. The image / is recovered by the application of 
{Fi,Fo)-^F-o to each dimension of /. Accurately recovering I from I requires the condition number 
of F^Fd to be small. The maximal condition number of Fd, which was obtained while applying 
Algorithm to various sizes n, is illustrated in Fig. The reason the term “maximal” is used, is 
since during the running of the algorithm, several operators Fd are used for a given n. As can be 
seen, the condition number is less than 25 even for very large volumes. The operator Fh is applied 
efficiently to a vector w by adding q — 1 zeros between every two samples of w, applying the inverse 
FFT to the resulting vector, and keeping the n central elements of the FFT-transformed vector 
(see Algorithm for a detailed description). The pseudo-code of the algorithm for recovering the 
volume from its samples in ffc is described in Algorithm 


Algorithm 3 AdjFDecimated: Applying F/, of Eq. 


2 T 1 to a vector 


1 : Input: y - Input vector of odd length. 
2 : g - integer (decimation factor). 

3: Output: X = F^y 
4: n <(—length(|/) — 1 
5: m ^ qn + 1 
6: z ^ zeros(m, 1) 

7: z{l : q : m) y 
8 : I length(^) 

9: pf <(— floor(//2) 

10 : Pc ^ ceil(//2) 

11: X ■(- z{[pf + 1:1, 1 : p/]) 

12 : xjfft <^IFFT(a:) 

13: a: ■(- m • Xifft{[Pc + 1:1, 1 : Pc]) 

14: X +- x{n + 1 : 2n) 
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Algorithm 4 InvDecimatedFreq. This procedure inverts the decimated frequencies and restores 
I from I of Eq. 


Eq. 
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Input: I - Matrix of size (n + 1) x (n + 1) x (n + 1). 


Output: I - Volume of size n x n x n such that Eq. 22 holds, 
c -ir- zeros(n, 1) 
m ^ qn + 1 

for k = —n/2 ,..., n/2 — 1 do 
for I = —u/2,... ,n/2 do 

c{k + n/2 + l) ^ c{k + n/2 + 1) + 

end for 
end for 

(Ml, M2, M3, M4) Toeplitzlnv(c) 

Di ToeplitzDiag(Mi), i = 1 ,..., 4 
Ji zeros(n + 1, n + 1, n) 

I2 ^ zeros(n + 1, n, n) 

I zeros(n, n, n) 
for fc = 1,..., n + 1 do 
for Z = 1,..., n + 1 do 

V I{k^ I ,:) 

V t—AdjFDecimated(u) 

u ^ToeplitzInvMul(T)i, D2, D^, D^, v) 

/i(fc, I,:) u 

end for 
end for 

for fc = 1,..., n + 1 do 

for / = 1,..., n do 

V Ii{k, :,l) 

V AdjFDecimated(u) 

u ToeplitzInvMul(T)i, D2, -D3, D4, v) 

-^2(^5 ^5 0 u 

end for 
end for 

for fc = 1,..., n do 
for / = 1,..., n do 

V ■(— /2(:, k, 1 ) 

V ■(—AdjFDecimated(u) 

u ^ToeplitzInvMul(Di, D2, D^, D4, v) 

I{:,k,l) t— M 

end for 
end for 


> Algorithm 6 

> Algorithm 7 


> Algorithm!^ 

> Algorithm^ 
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Condton nunt>er vs size 



Figure 6: The maximal condition number of Fp (see Eq. 21) obtained while applying Algorithm]^ 
to various sizes n. 


3.3 Complexity Analysis 

The inversion Algorithm consists of two steps: Resampling (lines 1-12) and recovering / from I 
on flc (line 13). The hrst step, resampling, which utilizes Toeplitz matrices (Algorithm]^, takes 
0{n\ogn) operations when preprocessing took place. This function (Toepliz-based resampling) is 
called 0{n) times in the recover2d function (Algorithm]^, which in turn is called also 0{n) times 
by Algorithmic A total of 0{n^\ogn) operations are needed for the resampling step. 

The second step, recovering I from I on 12^5 takes 0{'n?\ogn) operations. This yields a total 
computational complexity of 0{nHogn) operations. If the preprocessing is done in run time, 
the total complexity becomes 0{n'^) operations, due to the Durbin-Levinson algorithm used for 
inverting Toeplitz matrices. The computational complexity can be reduced from 0{n^) operations 


to 0{n^\og^n) by solving Eq. 11 as described in 

The total storage requirement for n applications of Algorithm is This storage is used to 
store the diagonal Toeplitz matrices Di,..., D/^. 


4 Numerical results 

We implemented Algorithm [C in MATLAB and applied it to several volnmes of sizes n x n x n 
where n = 32,64, ...,512. All the experiments were execnted on a Linux machine with two Intel 
Xeon processors (CPU X5560) running at 2.8GHz, with 8 cores in total and 96GB of RAM. All 
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Computation lir 



Figure 7: Runtime comparison between LS-based resampling and Toeplitz-NUFFT-based resam¬ 
pling (Algorithm for volumes of sizes n x n x n, consisting of PPFT applied to a normally 
distributed i.i.d. random samples with zero mean and unit variance. 


the 3D experiments were performed with q = 3 (see Eq. [^. 

Algorithm is based on a series of one-dimensional resampling operations. We compare two 
methods for implementing this one-dimensional resampling - least-squares-based (LS) approach 
and Toeplitz-NUFFT-based approach, both described in Section 2^ The LS-based approach, de¬ 
scribed in Eq. [T^ consists of Ending the coefficients vector a of a trigonometric polynomial, followed 
by direct evaluation of the polynomial at the resampling points. Using this resampling approach in 
Algorithm results in computational complexity of 0{n^) operations (excluding a preprocessing 
step). Nevertheless, its implementation in MATLAB is highly optimized. The Toeplitz-NUFFT- 
based approach, described in Section 2^ and Algorithmic results in computational complexity of 
Algorithm IC of (!7(?7 ,^logn-|-log 1/e) operations. The two resampling methods are compared by 
using them as the underlying one-dimensional resampling in Algorithm [C Both algorithms use a 
preprocessing step which is excluded from the reported running time. The error incurred by the 
two algorithms is measured by 




\ 




where I[k,m,n] is the original volume and Ir[k,m,n] is the reconstructed volume. Running times 
for both methods for volumes of sizes n x n x n, where n = 32,64, ...,512, are shown in Fig. 
The volumes in this experiment consist of PPFT transformed random independent samples from a 
normal distribution with zero mean and unit variance. Figure |C compares between the reconstruc¬ 
tion errors of Algorithm |C and those of LS-based resampling. The results show that both methods 
are comparable in terms of accuracy. 

Since it is currently impossible to process 3D volumes of sizes larger than n = 512, we 
compare between the LS-based resampling and the Toeplitz-based resampling by applying both 
one-dimensional resampling methods to a one-dimensional Chirp signal defined by cos(10A:|), 
kj = —TT -f i = 0, ...,n, n = 512, 1024, 2048, 4096. The signal’s samples given at kj 
were resampled to 0.3fcj using LS-based resampling and Toeplitz-based resampling. The one- 
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Toeplitz resampling vs least square resampling 



Figure 8: Error comparison between LS-based resampling and Toeplitz-NUFFT-based resampling 
(Algorithm]^ for volumes of sizes nxnxn, consisting of PPFT applied to a normally distributed 
random samples with zero mean and unit variance. 


dimensional Chirp signal is displayed in Fig. The running time of the two methods, which does 
not include pre-processing timing, appears in Fig. 10 The approximation errors of both methods 


are practically identical as shown in Fig. 11 


Next, we compare between the direct inverse PPFT in Algorithm [T| the iterative algorithm 
described in [35], and a single iteration of an implementation that computes the Gram operator 
of the PPFT using the NUFFT as suggested by [13]. For the latter method, several iterations are 
required, whose number depends on the condition number of the operator. However, to keep all 
timings within the same scale, we compare the other two algorithms to only a single iteration of 
the latter one. The results are illustrated in Fig. 1^, which shows that Algorithm [T] is faster than 
the iterative algorithm [39] as well as than a single iteration of the 3D NUFFT-based algorithm 
proposed by [T3]. Results for n = 512 do not appear in Fig. 12 (unlike Fig. as it was impossible 
to process volumes of that size using [39l fT3] . 

Given the long running time of the 3D NUFFT-based algorithm, we executed above only a single 
iteration of this algorithm. However, to compare accuracy of our proposed method with that of 
the 3D NUFFT-based algorithm, we execute next the 3D NUFFT-based algorithm until the error 
becomes smaller than 10“^^ or the number of iterations exceeds 100. Due to time constraints, we 
used only n = 16, 32, 64. The input for each case was of size n x n x n oi random normally 
distributed i.i.d. samples with zero mean and unit variance. We implemented the 3D NUFFT- 
based algorithm with the same preconditioner as in [39]. For n = 16 the algorithm took 112 
seconds and the resulting error was 7.25 x 10“^^, obtained after 40 iterations; for n = 32, it took 
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Figure 9: One-dimensional Chirp signal, used for testing the LS-based resampling and the Toeplitz- 
NUFFT-based resampling. 



Figure 10: Runtime comparison between the LS-based resampling and the Toeplitz-NUFFT-based 
resampling (Algorithm for a one-dimensional Chirp signal of length n. 
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Figure 11: Error comparison between the LS-based resampling and Toeplitz-NUFFT-based resam¬ 
pling (Algorithmic for a one-dimensional chirp signal of length n. 


Coinpuijfan time vs votume see 



Figure 12: Comparison between the running time of Algorithmic the iterative algorithm [SH], and 
the NUFFT-based algorithm [T3] . 
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(b) Bird, n = 128 
Figure 13: Test volumes. 



(c) Cup, n = 256 


Name 

Volume 

Time [sec] 

RMSE 

Dolphin 

64 X 64 X 64 

7.43 

1.69 X 10-^^ 

Bird 

128 X 128 X 128 

40.15 

3.6 X 10-^® 

Cup 

256 X 256 X 256 

270 

1.25 X 10-1^ 


Table 1: Results for volumes of Fig. 


13 


2, 710 seconds with error of 7.5 x 10“^^ after 100 iterations; for n = 64 it took 22,477 seconds (more 
than 6 hours) with error of 1.61 x 10“^^ after 100 iterations. The iterative algorithm described 
in [39] gives similar errors to the direct inversion algorithm (Algorithm whose error appears in 

Fig. ^ 

Finally, we tested the algorithm on real volumes of different sizes: a dolphin of size 64 x 64 x 64 


(Fig. 13a), a bird of size 128 x 128 x 128 (Fig. 13b), and a 3D cup of size 256 x 256 x 256 (Fig. 13c). 
The volumes were taken from the McGill 3D shape dataset and PPFT was applied to each 
one them to be used as an input for the inversion algorithm. The results appear in Table 


5 Conclusions 

In this paper, a new algorithm for inverting the 3D pseudo-polar Fourier transform (PPFT) is 
described. The algorithm processes at each iteration a two-dimensional slice of the input, where 
each such processing uses only one-dimensional operations. The main component of the algorithm 
is fast resampling of univariate trigonometric polynomials. The resampling is implemented using 
a ID non-uniform Fourier transform together with fast algorithms for Toeplitz matrices. The 
algorithm is not iterative, and requires a fixed amount of time that depends only on the size of 
the input. Moreover, the algorithm has low memory requirements, allowing to process large 3D 
datasets in a reasonable time. The performance of the algorithm is demonstrated on volumes as 
large as 512 x 512 x 512 in double precision. 
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Appendices 

A Toeplitz solvers algorithms 


Algorithm 5 ToeplitzResample; Fast resampling of a univariate trigonometric polynomial. 
1: Input: y - vector of length n. 

2: / - values of the trigonometric polynomial / on the set y. 

3: X - vector of length n. 

4: Output: g - values of / on the set x. 

5 : fc [—n /2 : n /2 — 1 ] 

6: c length(|/) x NUFFTi(|/, 

7 : [Ml, M2., Mi) <(— ToeplitzInv(c, c) 

8 : for i = 1 ,..., 4 do 

9 : Di ^ToeplitzDiag(Mi(:, 1 ), Mj(l,:)) 

10: end for 

11: V •(— length {y)x NUFFTi(?/, /, —1, n) 

12: a •(—ToeplitzInvMul(Zi)i, D2, D^, D 4 , v) 

13: g NUFFT2 {x, 1, n, a) 


> fc is a row vector 
> k{l) is the first element of k 

> Can be computed in the preprocessing step 

> Can be computed in the preprocessing step 


NUFFTfc (fc = 1, 2) refers to the type of the NUFFT - see Eqs. 15 and 16 The parameters —1 
and 1 in lines 6,11,13 in Algorithm refer to the sign of i in the complex exponent. The function 
Toeplitzinv is described in Algorithm and computes the Gohberg-Semencul factorization of a 
symmetric Toeplitz matrix. 


Algorithm 6 Toeplitzinv: Factorize a symmetric Toeplitz matrix using the Gohbcerg-Semencul 
formula. _ 

1: Input: c - First column (row) of a symmetric Toeplitz matrix. 

2: Output: Ml, M2, M3, M4 - Gohberg-Semencul factorization. 

3: n <(—length(c) 

4: Co ^ zeros(n, 1) 

5: Cn ^ zeros(n, 1) 

6: eo(l) ^ 1 
7: e„(n) ^ 1 
8: Solve T„(c, c)x = cq 
9 : Solve T„(c,c)i/= Cn 

10: Gompute Mi using Eq. 

11: Gompute M 2 using Eq. 

12: Gompute M3 using Eq. 

13: Gompute M4 using Eq. 
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Algorithm 7 ToeplitzDiag; Diagonalize the circular embedding of a Toeplitz matrix 

1 : Input: c, r - First column/row of the Toeplitz matrix. 

2 : Output: D - Diagonal form of Tn{c,r) embedded in a circulant matrix. 

3: n <(— length(c) 

4: C ^ [c; 0; r {n : —1 ; 2)] 

5: D ^ FFT(C') 


Algorithm 8 ToeplitzInvMul: Mulitply an inverse Toeplitz matrix A ^ of size n x n by a vector 
V, given the diagonal forms of the Gohberg-Semencul factors of A~^ in 0{n\ogn) operations. 

Input: Di, D2, D^, D4 - Diagonal forms of the Gohberg-Semencul factors of the inverse Toeplitz 
matrix. 

V - vector of length n. 

Output: u - The product A~^v. 

Ui ^ToeplitzMul(D2, n) 

Ml ^ToeplitzMul(Di,Ml) 

M2 ^ToeplitzMul(D 4 , m) 

M2 ^ToeplitzMul(D3, M2) 

M ■(— Ml — M2 


Algorithm 9 ToeplitzMul: Multiply a Toeplitz matrix A of size n x n, whose diagonal factor is 
D, by a vector v in 0{nlogn) operations. 

1 : Input: D - Diagonal factor of A (computed using Algorithm]^. 

2: V - vector of length n. 

3: Output: M - the product Av. 

4: n length (m) 

5: M •(— [m; zeros(n, 1)] 

6 : m; FFT(m) 

7: u IFFT(D. * w) 

8: M m(1 : n) 
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